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We update our previous work of a Langevin-with-interaction model for charmonium in heavy- 
ion collisions, by considering the effect due to recombination. We determine the contribution to 
J/ip yields from cc pairs whose constituent quarks originate from two different hard processes. 
Like the surviving J/ip states, the recombinant J/ip also undergo both a stochastic interaction, 
determined by a hydrodynamical simulation of the heavy-ion collision, and an interaction determined 
by the QQ potentials measured on the lattice for appropriate temperatures. From the results 
of these simulations, we determine both the direct and the recombinant contribution to the J/ip 
yields for RHIC conditions, and find that for central collisions, between 30% and 50% of the J/ip 
yield is due to recombinant production. We compare our results with other models and look for 
how the recombinant contribution differs from the surviving contribution in the differential pr 
yields. Including the recombinant contribution improves the agreement with the latest analysis of 
charmonium at RHIC, which shows an absence of anomalous suppression except in the most central 
collisions. 



PACS numbers: 

I. INTRODUCTION 

In a previous paper [1] , we argued that the microscopic 
dynamics of charmonium in a heavy-ion collision should 
be modeled as a stochastic interaction with strongly- 
coupled quark-gluon plasma (sQGP). In such a plasma 
the diffusion coefficient is very small, leading to rapid 
thermalization in momentum space and slow spatial dif- 
fusion, which is further slowed by the attraction of the 
constituent charm and anti-charm quarks in the pair. We 
concluded that the amount of equilibration in space (and 
thus the J/ip yield) depends strongly on the timescales of 
the collision. In realistic simulations of Au+Au collisions 
at RHIC, where the sQGP phase exists for r ~ 5 fm/c , 
our model predicted a survival probability for J/ip ~ 1/2 
even in the most central collisions, much larger than pre- 
viously expected. Those results were able to explain qual- 
itatively the data from PHENIX. 

In pQ, we did not consider another possibly important 
source of charmonium at RHIC, the "recombinant" con- 
tribution of J/ip particles whose constituent quarks origi- 
nate from different hard processes. At very high collision 
energies, as the number of charm pairs per event grows, 
recombinant charmonia could potentially lead to an en- 
hancement of the final J/ip yields in a heavy-ion collision, 
reversing the current suppression trend. Using the grand 
canonical ensemble approach, Braun-Munzinger and col- 
laborators m determine the fugacity of charm by the 
number of cc pairs produced initially. The "statistical 
hadronization" approach to charmonium assumes com- 
plete thermal equilibration of charmonium. Another ap- 
proach has been taken by Grandchamp and Rapp |4J , who 
treat the J/ip yields from heavy-ion collisions as coming 
from two sources: the direct component, which decays 
exponentially with some lifetime Ta, and the coalescent 
component, which is determined by the same mechanism 



in [2], with the additional consideration that spatial equi- 
libration of charm does not happen. To account for en- 
hanced local charm density, by small spatial diffusion, 
they had introduced another factor - the "correlation vol- 
ume" V corr - which was estimated. The present work can 
be viewd as a quantitative dynamical calculation of this 
parameter. 

To gain insight, we should compare these models with 
our model in [1]. The Langevin-with-interaction model 
for cc pairs in medium makes no assumptions about 
complete thermalization, and shows how even in central 
Au+Au collisions at the RHIC, the majority of the J/ip 
yields may survive the QGP phase. However, the model 
predicts rapid thermalization in the momentum distribu- 
tions of charmonium, as well as equilibration in the rel- 
ative yields of the various charmonium states due to the 
formation of "quasi-equilibrium" in phase space. This 
requires no fine-tuning of the rates for charmonium in 
plasma; it is just a natural consequence of the strongly 
coupled nature of the media, detailed by the Langevin 
dynamics. However, recombinant production of charmo- 
nium may still be an important effect in our model, due 
to the fact that in central collisions, the densities of un- 
bound charm quarks can be quite high in some regions 
of the transverse plane. 

Our model simulates an ensemble of cc pairs, gen- 
erated initially by PYTHIA event generation and 
then evolved according to the Langevin-with-interaction 
model. We evolve the pairs not assuming any form of 
equilibrium, and then average over possible pairings of 
the quarks to form recombinant charmonium. 

The outline of this work is as follows: in Section HT1 we 
will describe how we simulated charm in plasma and took 
into account the contribution due to recombinant J/ip , 
and in Section |lTT| we take the opportunity to describe the 
progress in this model so far and also to summarize where 
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future work is needed for a Langevin-with-interaction de- 
scription of J/ip suppression. In Appendix |A| we discuss 
the statistics necessary to calculate the recombinant con- 
tribution to J ftp yields. 



II. RECOMBINANT CHARMONIUM IN 
HEAVY-ION COLLISIONS 

A. Langevin-with-interaction model for cc pairs in 
a heavy-ion collision 

As we have done in our previous paper, we simulate 
J ftp in medium with a hydrodynamical simulation of the 
collision. As before, we start with a large ensemble of 
cc pairs whose momenta are determined with PYTHIA 
event generation [3 : . The positions of the initial hard col- 
lisions in the transverse plane at mid-rapidity are deter- 
mined by sampling the distribution in N co u determined 
from the Glauber model. In this way, our local densities 
of cc pairs vary as one would expect from the Glauber 
model, which gives an enhancement for recombination 
towards the center of the transverse plane. Each element 
of the ensemble now contains N cc pairs. The number of 
pairs N depends on the impact parameter of the collision 
and needs to be determined. 

The average number of cc pairs for a Au+Au collision 
at RHIC varies with impact parameter and has been in- 
vestigated by the PHENIX collaboration at mid-rapidity 
[B]. The measured cross sections for charm production 
vary somewhat with the centrality of the collision and 
achieves a maximum of about 800 fib for semi-central 
collisions. The nuclear overlap function Taa{V) can be 
calculated with the Glauber model. We used a conve- 
nient program by Dariusz Miskowiec \J\ to evaluate this 
function. With a centrality dependent cross-section cj Ec , 
we can easily calculate the average number of cc pairs 
in a collision: N Sc = Taa&cc- The number of cc pairs 
reaches a maximum in central collisions, with an average 
of 19 pairs per collision. 

In order to determine the probability for two charm 
quarks from different hard processes to form recombinant 
charmonium, we must average over the different possible 
pairings of all of the unbound quarks in each element of 
our ensemble. This is discussed in Appendix |A"| in gen- 
erality Since the number of cc pairs approaches 20 for 
central Au+Au collisions at RHIC, we are faced with an- 
other issue: there are 20! possible pairings and it has be- 
come impractical to calculate the probability of each in- 
dividual pairing this way. In general, we would be forced 
to perform permutation sampling of this partition func- 
tion, preferably with some Metropolis algorithm. How to 
sample over permutations with a Metropolis algorithm is 
discussed thoroughly in the literature, for an excellent 
review of this see Ceperley j5j. However, for RHIC, the 
situation simplifies due to the relatively low densities of 
cc pairs involved. We ran our simulation for the most 
central Au+Au collisions at RHIC and examined how 



many "neighbors" each charm quark had. A "neighbor" 
is defined as a charm anti-quark, originating from a dif- 
ferent pQCD event yielding the given charm quark, which 
is close enough to the charm quark that it could poten- 
tially form a bound state, in other words r is such that 
V corne u(r) < 0.88 GeV. The number of charm quarks 
expected to have one and only one neighbor in the most 
central Au+Au collisions was found to be 5.5%, while 
only 0.2% of the charm quarks are expected to have more 
than one neighbor. Therefore, even at the most central 
collisions at RHIC, we can be spared possibly compli- 
cated permutation samplings. Of course, this situation 
is not true in general, and for the numbers of pairs pro- 
duced in a typical heavy-ion collision at the LHC one 
should modify these combinatorics. 



B. New analysis of the data including improved 
dAu sample 

The data with which we now compare our results is 
different from that which we used for comparison in our 
previous work. New data analysis of Au+Au and d+Au 
described in [10] account for the (anti-)shadowing and 
the breakup of charmonium due to the cold nuclear mat- 
ter effects (parameterized by <7 a bs) hi the framework of a 
Glauber model for the collision. The calculations at for- 
ward and mid-rapidity are now done independently, since 
shadowing and breakup could be considerably different at 
different rapidities. This new analysis is a significant suc- 
cess, demonstrating the high suppression at forward ra- 
pidity (previously very puzzling) as being due to cold nu- 
clear matter effects. New ratios of observed suppression 
due to cold nuclear matter Raa/ Raa M > plotted versus 
the energy density times time er, show common trends 
for RHIC and SPS data, which was not the case previ- 
ously. We use this new analysis as a measure of survival 
probability in our calculation. 



C. The results 

Before we show the results, let us remind the reader 
that our calculation is intended to be a dynamical one, 
with no free parameters. We use a hydrodynamical sim- 
ulation developed in [8 which is known to describe accu- 
rately the radial and elliptic collective flows observed in 
heavy-ion collisions. Our drag and random force terms 
for the Langevin dynamics has one input - the diffu- 
sion coefficient for charm - constrained by two indepen- 
dent measurements (pt distributions and V2(pt) mea- 
surements for single lepton - charm - performed in Ref. 
[S] . The interaction of these charm quarks are determined 
by the correlators for two Polyakov lines in 1QCD [11 . 

Having said that, we still are aware of certain uncer- 
tainties in all the input parameters, which affect the 
results. In order to show by how much the results 
change if we vary some of them, we have used the un- 
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FIG. 2: (Color online.) The surviving and recombinant J/ip 
yields, plotted versus the radial distance from the center of 
the transverse plane. 



stroyed by the medium, which ultimately coalesce in the 
end. Second, the direct J/ip states' abundance, when 
compared with the abundances of excited charmonium 
states, does not necessarily need to be as expected from 
these particles' Boltzmann factors. For our model, these 
relative abundances do make sense for direct charmonium 
states, due to the formation of a quasi-equilibrium dis- 
tribution. 



FIG. 1: (Color online.) R A ^ malous = Raa/RaT' for J/ip 
versus centrality of the AuAu collisions at RHIC. The data 
points with error bars show the PHENIX Au+Au measure- 
ments with cold nuclear matter effects factored out as in 
[10] . Other points, connected by lines, are our calculations 
for the two values of the QCD phase transition temperature 
T c = 165 MeV (upper) and T c = 190 MeV (lower). From bot- 
tom to top: the (green) filled squares show our new results, 
the recombinant J/ip , the open (red) squares show the Raa 
for surviving diagonal J/ip , the open (blue) circles show the 
total. 



certainty in the value for the critical temperature T c . 
For these reasons, we show the results for two values 
T c = 165, 190 MeV, in FigJTJ 

As can be seen, a higher T c value improves the agree- 
ment of our simulation with the latest analysis of the 
data, because in this case the QGP phase is shorter in 
duration and the survival probablity is larger. However 
the recombinant contribution (shown by filled squares) 
is in this case relatively smaller, making less than 1/3 of 
the yield even in the most central collisions at RHIC. 

Our results for the total, direct, and recombinant con- 
tributions resembles considerably the results of Zhao and 
Rapp obtained from their two-component model |12j . 
However it is important to point out two important differ- 
ences. First of all, what is described by Zhao and Rapp 
as the second component due to statistical coalescence 
includes, with the recombinant J/ip , surviving J/ip , de- 



D. Recombinant J/ip and ^-distributions 

So far, we have only considered the effect of the recom- 
binant production on the overall yields of J/ip particles 
at the RHIC. We should test our model by considering 
whether or not adding the recombinant contribution can 
change the shape of any differential J/ip yields. 

One differential yield where we may expect the surviv- 
ing and recombinant component to have different behav- 
iors is in the p^-distributions for central Au+Au colli- 
sions. The surviving J/ip states tend to originate in the 
periphery of the collision region, since the J/ip states pro- 
duced here endure the sQGP phase for the shortest times. 
However, the recombinant contribution should form to- 
ward the center of the collision region, since this is where 
the density of initial cc pairs is highest, and as we have 
been showing for some time, spatial diffusion is incom- 
plete in the sQGP. Therefore, since the effect due to flow 
on the pr-distributions has Hubble- like behavior, with 
the radial velocity of the medium scaling with distance 
from the center of the transverse plane, we would ex- 
pect the recombinant contribution to exist, on average, 
in regions of the medium with significantly smaller flow 
velocities. 

Figure [2] demonstrates this behavior existing in our 
simulation. 

We should now determine whether or not this differ- 
ence of the yield versus r can be observed in the J/ip 
yield versus pt- As we have shown in our previous paper, 
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FIG. 3: (Color online.) The surviving and recombinant J/ip 
yields versus pt- 



during the phase transition from QGP to the hadronic 
phase in heavy- ion collisions, our model predicts a small 
change in the total J/i/j yield but relatively large changes 
in the J/ip pt distributions, with these changes strongly 
dependent on the drag coefficient for quarkonium during 
this time, and T c . We can easily run our code with an 
LH8 equation of state and make several predictions for 
the two components' p t distributions. However, for rea- 
sons which will become apparent, we are only interested 
in the upper limit of the effect of flow on p t in Au+Au 
collisions at the RHIC. Therefore, we ran our simulation 
where we assumed a phase transition which lasts 5 fm/c, 
during which the J/ip particles have a mean free path of 
zero, in a Hubblc-likc expansion. 

The p t distributions after this expansion are shown in 
Figure [3j It is visible from this plot that the recombinant 
contribution will observably increase the total yield at 
low p t (where the total yield is significantly higher than 
the surviving component alone) and have little effect at 
higher p t (where the total and the surviving component 
alone are nearly the same). However, we have found that 
even in this extreme limit, there is no clear signal in the 
differential pt yields for there being two components for 
J/t/j production at the RHIC. 

This test, however, should not be abandoned for mea- 
surements of the differential yields at higher collision en- 
ergies. Since the recombinant contribution grows sub- 
stantially as charm densities are increased, it should be 
checked whether or not the recombinant contribution is 
more strongly peaked in the center of the transverse plane 
of LHC collisions, and whether or not two components 



to the differential yields should become observable there. 
We will follow up on this issue in a work we have in 
progress. 



III. DISCUSSION 

We have found that at central Au+Au collisions at 
RHIC the fraction of recombinant pairs should be con- 
siderable, up to 30-50%, with smaller fractions at more 
peripheral collisions. The exact number depends on de- 
tails of the model, such as duration of the QGP phase 
and the magnitude of the critical temperature T c . We 
have also gone a step further, and attempted to find dif- 
ferent behaviors of these two components in differential 
yields, so that these two components might be disentan- 
gled. This test (examining the differential p t yields) fails 
to identify clearly two different components. We will pur- 
sue whether or not this test works for the yields at the 
LHC. 

Our model for charmonium in sQGP is rather con- 
servative: we merely assume that the constituent charm 
quarks experience dynamics similar to the Langevin dy- 
namics of single charm quarks in SQGP, which has al- 
ready shown good agreement with the Raa{pi) and 
v-2{pt) measured at PHENIX for single charm. 

One final, careful observation of our results is worth 
mentioning. As one can see from our results of Fig. [l] the 
model seems to be working well for central collisions, in 
which there is a QGP phase lasting for several fm/c and 
leading to a possibility for charm quarks to diffuse away 
from each other, far enough so that J/ip states would 
not survive. However it overpredicts suppression for pe- 
ripheral collisions, which if the cold matter analysis 
will hold against further scutiny - is nearly completely 
absent. One possible reason for that can be survival of 
the flux tubes (QCD strings) between quarks well into 
the mixed phase or even in small region of temperatures 
above T c , as was recently advocated by one of us [13 in 
connection with "ridge" phenomenon. 
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APPENDIX A: CANONICAL ENSEMBLES FOR 
N cc -PAIR SYSTEMS 

In this section we will determine a partition function 
for a canonical ensemble of N charm pair systems (that 
is, an ensemble of very many systems, where each system 
contains N cc pairs) which correctly averages over differ- 
ent possible pairings of charm and anticharm quarks and 
can therefore describe recombination in heavy-ion colli- 
sions. This averaging is possible computationally but is 
non-trivial, and for RHIC collisions we will take a binary 
approximation which makes this averaging much easier. 
We argue, however, that the unsimplified approach is 
necessary for describing collisions at the Large Hadron 
Collider, and for this reason we include this discussion 
here. 

Our simulation could be thought of as a canonical en- 
semble description of charmonium in plasma: we can 
think of our large set of cc pairs as a set of systems, each 
system containing N pairs, with each system's dynam- 
ics modeled as a stochastic interaction with the medium, 
with a deterministic interaction of each heavy quark with 
the other quark in the pair. Each system in this set 
samples the distribution of cc pairs in the initial col- 
lision, the geometry of the collision, and also samples 
the stochastic forces on the heavy quarks. Up to this 
point, we have only thought of each system of this set as 
consisting of a single cc pair. The interaction of charm 
quarks from different hard events is negligible compared 
with the stochastic interaction and the interaction within 
the pair, partly because near T c , the dynamics of charm 
pairs seems best described with some generalization of 
the Lund string model, which allows no interaction be- 
tween unpaired charm quarks [3], Therefore, it is simple 
bookkeeping to think now of the systems as each consist- 
ing of N cc pairs. 



However, even though the dynamics of the system is 
not changed when considering many cc pairs per colli- 
sion, the hadronization ("pairing") of these 2N charm 
quarks is now a non-trivial issue. For simplicity, assume 
that the quarks all reach the freezeout temperature T c 
at the same proper time. There are Nl different possible 
pairings of the quarks and anti-quarks into charmonium 
states (each pairing is an element of the permutation 
group Sn)- Call a given pairing a (which is an element 
of Sjv)- Near T c , the relative energetics of a pairing a is 
given by 



E, 
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where V(r) is the zero-temperature Cornell potential 
(with some maximum value at large r, corresponding to 
the string splitting), Fj the position of the i-th charm 
quark, r the position of the z-th charm antiquark, and 
a(i) the integer in the i-th term of the permutation. 

One way to proceed is to average over these pairings 
according to their Boltzmann factors. In this way, the 
probability of a given pairing would be given by 



p(i) 



i m 

- exp(-Ei/T c ), Z = £exp(-£ 4 /T c ) 



(A2) 



However, this averaging ignores the possibility of sur- 
vival of bound cc states from the start to the finish of 
the simulation, in that pairings which artificially "break 
up" bound states are included in the average. This goes 
against the main point of our last paper: that it is ac- 
tually the incomplete thermalization of cc pairs which 
explains the survival of charmonium states. 

For this reason, the averaging we perform rejects per- 
mutations which break up pairs that would otherwise be 
bound: we average over a subgroup S N of Sn, and de- 
termine the probability based on this modified partition 
function: 

P(i) = ^exp(-E t /T c ) 1 Z=Y1 exp(-^/T c ), (A3) 

<?es' N 

where Ei specifies the energy of a pairing we permit. We 
will average over the permutations in this way. 

By doing this, we will use a fully canonical ensem- 
ble description for charm in plasma, which holds for any 
value for N, large or small. Previous work in statistical 
hadronization used the grand canonical approach to ex- 
plain relative abundances of open and hidden charm [5] , 
which can only be applied where thermalization may be 
assumed to be complete. 



